# packages

library(tidyverse)
library(margins)
library(marginaleffects)
library(pROC)
library(caret)
library(ResourceSelection)
library(generalhoslem)
library(modelsummary)
library(stargazer)
library(countrycode)
library(peacesciencer)
library(WDI)
library(tidytable)
library(broom.mixed)
library(lmtest)
library(sandwich)
library(effects)
library(Zelig)


# General Data Preparation


pmcdata<-Petersohn_et_al_CMA_Dataset_Version2
pmcdata<-pmcdata%>%filter(pmcdata$Clnt==1|pmcdata$Clnt==2|pmcdata$Clnt==5)
summary(pmcdata)
pmcdata[pmcdata==999]<-NA
summary(pmcdata)
pmcdata[pmcdata==998]<-NA
summary(pmcdata)
table(pmcdata$CMA)
table(pmcdata$Serv)
attach(pmcdata)
pmcdata$clientcode<-countrycode(pmcdata$ClNa,origin = "cowc",destination = "cown")
table(pmcdata$COWCCode)
table(pmcdata$clientcode)
length(pmcdata$COWCCode)
fliter1<-unique(pmcdata$COWCCode)
fliter2<-unique(pmcdata$clientcode)
unique(pmcdata$clientcode)


# Getting Data From ATOP and COW and World Bank


# COW and ATOP data


cow_data<-create_dyadyears(system = "cow",mry = T,subset_years = c(1980:2016),directed = F)%>%add_nmc()%>%add_atop_alliance()
cow_data1<-cow_data%>%dplyr::select(c(ccode1,upop1,tpop1,irst1,cinc1,pec1,milper1,milex1,year,atop_defense))
cow_data2<-cow_data%>%dplyr::select(c(ccode1,ccode2,upop1,tpop1,irst1,cinc1,pec1,milper1,milex1,year,atop_defense))
cow_data2$nodefwithclnt<-ifelse(cow_data2$ccode1==fliter1&cow_data2$ccode2!=fliter2,1,0)
cow_data2$defwithclnt<-ifelse(cow_data2$ccode1==fliter1&cow_data2$ccode2==fliter2,1,0)
table(cow_data2$defwithclnt)
cow_nmc11<-cow_data2%>%dplyr::filter(ccode1==fliter1)
cow_nmc1<-cow_data1%>%dplyr::filter(ccode1==fliter1)
names(cow_nmc1)[9]<-"YEAR"
names(cow_nmc11)[10]<-"YEAR"
cow_nmc2<-cow_data1%>%dplyr::filter(ccode1==fliter2)
cow_nmc22<-cow_data2%>%dplyr::filter(ccode1==fliter2)
names(cow_nmc2)[9]<-"YEAR"
names(cow_nmc22)[10]<-"YEAR"
cow_nmc2<-subset(cow_nmc2,select = -atop_defense)
cow_nmc22<-subset(cow_nmc22,select = -c(atop_defense,nodefwithclnt,defwithclnt))
newdata<-left_join(pmcdata,cow_nmc1)
newdata1<-left_join(pmcdata,cow_nmc11)
newdata2<-distinct(newdata,EchangeID,.keep_all = T)
newdata21<-distinct(newdata1,EchangeID,.keep_all = T)
newdata2<-subset(newdata2,select = -ccode1)
newdata21<-subset(newdata21,select = -c(ccode1,ccode2))
summary(newdata2)
summary(newdata)
names(newdata2)[31]<-"urbpoptrgt"
names(newdata2)[32]<-"totpoptrgt"
names(newdata2)[33]<-"irontrgt"
names(newdata2)[34]<-"cinctrgt"
names(newdata2)[35]<-"econstrgt"
names(newdata2)[36]<-"milpertrgt"
names(newdata2)[37]<-"milextrgt"
names(newdata21)[31]<-"urbpoptrgt"
names(newdata21)[32]<-"totpoptrgt"
names(newdata21)[33]<-"irontrgt"
names(newdata21)[34]<-"cinctrgt"
names(newdata21)[35]<-"econstrgt"
names(newdata21)[36]<-"milpertrgt"
names(newdata21)[37]<-"milextrgt"
newdata3<-left_join(newdata2,cow_nmc2)
newdata31<-left_join(newdata21,cow_nmc22)
newdata3<-distinct(newdata3,EchangeID,.keep_all = T)
newdata31<-distinct(newdata31,EchangeID,.keep_all = T)
newdata3<-subset(newdata3,select = -ccode1)
newdata31<-subset(newdata31,select = -c(ccode1,ccode2))
table(newdata31$nodefwithclnt)
table(newdata31$atop_defense)
table(newdata31$defwithclnt)
summary(newdata31$defwithclnt)
names(newdata3)[39]<-"urbpopcln"
names(newdata3)[40]<-"totpopcln"
names(newdata3)[41]<-"ironcln"
names(newdata3)[42]<-"cinccln"
names(newdata3)[43]<-"econscln"
names(newdata3)[44]<-"milpercln"
names(newdata3)[45]<-"milexcln"
names(newdata31)[41]<-"urbpopcln"
names(newdata31)[42]<-"totpopcln"
names(newdata31)[43]<-"ironcln"
names(newdata31)[44]<-"cinccln"
names(newdata31)[45]<-"econscln"
names(newdata31)[46]<-"milpercln"
names(newdata31)[47]<-"milexcln"


# Mountain range of target


data("rugged")
rugged1<-rugged%>%filter(ccode%in%fliter1)
rugged1$newlmtnest<-exp(rugged1$newlmtnest)
rugged1<-subset(rugged1,select = -c(gwcode,rugged))
names(rugged1)[1]<-"COWCCode"
newdata3<-left_join(newdata3,rugged1)
newdata31<-left_join(newdata31,rugged1)
names(newdata3)[46]<-"mounttrgt"
names(newdata31)[48]<-"mounttrgt"


# Distance between Target and Client


distance<-create_dyadyears(system = "cow",mry = T,directed = F,subset_years = c(1980:2016))%>%add_capital_distance()
distance<-distance%>%filter(ccode1%in%fliter2&ccode2%in%fliter1)
names(distance)[1]<-"clientcode"
names(distance)[2]<-"COWCCode"
names(distance)[3]<-"YEAR"
maindata4<-maindata3
maindata4<-left_join(maindata4,distance)
maindata3<-left_join(maindata3,distance)


# GDP World bank


options(scipen = 999)
filter2.1<-na.omit(fliter2)
filter3<-countrycode(filter2.1,origin = "cown",destination = "iso3c")
gdp_client<-WDI(country = filter3,indicator = "NY.GDP.MKTP.CD",start = 1980,end = 2016)
names(gdp_client)[4]<-"YEAR"
names(gdp_client)[3]<-"GDP"
gdp_client$COWCCode<-countrycode(gdp_client$iso2c,origin = "iso2c",destination = "cown")
gdp_sub<-subset(gdp_client,select = c(GDP,YEAR,COWCCode))
names(gdp_sub)[3]<-"clientcode"
newdata3<-left_join(newdata3,gdp_sub)
newdata31<-left_join(newdata31,gdp_sub)
fliter1.1<-na.omit(fliter1)
filter4<-countrycode(fliter1.1,origin = "cown",destination = "iso3c")
gdp_trgt<-WDI(country = filter4,indicator = "NY.GDP.MKTP.CD",start = 1980,end = 2016)
names(gdp_trgt)[4]<-"YEAR"
names(gdp_trgt)[3]<-"GDPT"
names(gdp_trgt)[5]<-"COWCCode"
gdp_trgt$COWCCODE<-countrycode(gdp_trgt$iso2c,origin = "iso2c",destination = "cown")
gdp_sub<-subset(gdp_trgt,select = c(GDPT,YEAR,COWCCode))
maindata4<-maindata3
maindata4<-left_join(maindata4,gdp_sub)
maindata3<-left_join(maindata3,gdp_sub)


# War Data Set Interstate UCDP


cow_ucdp<-ucdp_acd
cow_ucdp<-cow_ucdp%>%dplyr::filter(gwno_a==fliter1)
cow_ucdp<-cow_ucdp%>%dplyr::filter(year>=1980&year<=2016)
cow_ucdp<-subset(cow_ucdp,select = c(incompatibility,year,intensity_level,type_of_conflict,ep_end,gwno_a))
names(cow_ucdp)[2]<-"YEAR"
names(cow_ucdp)[6]<-"COWCCode"
newdata3<-left_join(newdata3,cow_ucdp)
newdata31<-left_join(newdata31,cow_ucdp)


# Polity 2 Scores of Client and Target


polityscores<-p5v2018
polityscores<-polityscores%>%filter(year>=1980&year<=2016)
politytarget<-polityscores%>%filter(ccode%in%fliter1)
politytarget<-subset(politytarget,select = c(ccode,year,polity2))
names(politytarget)[1]<-"COWCCode"
names(politytarget)[2]<-"YEAR"
names(politytarget)[3]<-"polscrtrgt"
newdata3<-left_join(newdata3,politytarget)
newdata31<-left_join(newdata31,politytarget)
polityclient<-polityscores%>%filter(ccode%in%fliter2)
polityclient<-subset(polityclient,select = c(ccode,year,polity2))
names(polityclient)[1]<-"clientcode"
names(polityclient)[2]<-"YEAR"
names(polityclient)[3]<-"polscrclnt"
newdata3<-left_join(newdata3,polityclient)
newdata31<-left_join(newdata31,polityclient)


# Saving Data


library(xlsx)
write.xlsx(newdata3,"maindata1.xlsx",col.names = T,row.names = F,append = F)
write.xlsx(newdata31,"maindata2.xlsx",col.names = T,row.names = F,append = F)
write.xlsx(maindata4,"finaldata.xlsx",col.names = T,row.names = F,append = F)


# More exploration


table(Petersohn_et_al_CMA_Dataset_Version2$Serv)
table(maindata1$Serv)
xtabs(~maindata1$Serv+maindata1$Region,data = maindata1)
table(Petersohn_et_al_CMA_Dataset_Version2$Clnt)
table(maindata1$ForT)
xtabs(~ForT+Clnt,data = maindata1)
table(maindata1$Clnt)
table(maindata1$ForT)
summary(maindata1$ForT)
xtabs(~ForT+Cons+Clnt,data = maindata1)
xtabs(~ForT+Serv,data = maindata1)
cor(maindata1$ForT,maindata1$polscrclnt,use = "complete.obs")
glm1<-glm(as.factor(ForT)~polscrclt,data = maindata1,family = binomial(link = "logit"))
summary(glm1)
margins1<-margins::margins(glm1,type="response")
summary(margins1)


# Main Work on Data Set


summary(maindata2)
table(maindata2$ForT)
stargazer(as.data.frame(maindata2),type = "text",out = "Descriptive.txt",style = "ajps",align = T,digits = 5)
hist(maindata2$Serv)
table(maindata2$Serv)
xtabs(~ForT+Serv,data = maindata2)


# Modifying Main Data


table(as.numeric(as.factor(maindata2$incompatibility)))
maindata2$government<-ifelse(maindata2$incompatibility=="government",1,0)
maindata2$territory<-ifelse(maindata2$incompatibility=="territory",1,0)
table(maindata2$government)
table(maindata2$territory)
table(as.numeric(as.factor(maindata2$type_of_conflict)))
which(maindata2$type_of_conflict=="intrastate")
maindata2$intintra<-ifelse(maindata2$type_of_conflict=="II",1,0)
table(maindata2$intintra)
table(maindata2$intensity_level)
maindata2$forclnt<-ifelse(maindata2$ForT==0,1,0)
maindata2$forthird<-ifelse(maindata2$ForT==1,1,0)
table(maindata2$forclnt)
table(maindata2$forthird)
maindata2$combat<-ifelse(maindata2$Serv==1,1,0)
maindata2$mconsult<-ifelse(maindata2$Serv==2,1,0)
maindata2$mtrain<-ifelse(maindata2$Serv==3,1,0)
maindata2$mopsupport<-ifelse(maindata2$Serv==4,1,0)
maindata2$mlogist<-ifelse(maindata2$Serv==5,1,0)
maindata2$intellig<-ifelse(maindata2$Serv==6,1,0)
maindata2$secur<-ifelse(maindata2$Serv==7,1,0)
maindata2$sconsult<-ifelse(maindata2$Serv==8,1,0)
maindata2$strain<-ifelse(maindata2$Serv==9,1,0)
maindata2$slogist<-ifelse(maindata2$Serv==10,1,0)
maindata2$recons<-ifelse(maindata2$Serv==11,1,0)
table(maindata2$nodefwithclnt)
maindata2$defwithother<-ifelse(maindata2$atop_defense==1&maindata2$nodefwithclnt==1,1,0)
table(maindata2$defwithother)
table(maindata2$atop_defense)
table(maindata2$AgSt)
maindata2$groupPMSC<-ifelse(maindata2$AgSt==1,1,0)
maindata2$locPMSC<-ifelse(maindata2$AgSt==2,1,0)
maindata2$intPMSC<-ifelse(maindata2$AgSt==3,1,0)
maindata2$ambgPMSC<-ifelse(maindata2$AgSt==4,1,0)
maindata2$nontraded<-ifelse(maindata2$OwSt==1,1,0)
maindata2$traded<-ifelse(maindata2$OwSt==2,1,0)
maindata2$abspol<-abs(maindata2$polscrclnt-maindata2$polscrtrgt)
maindata2$logGDP<-log(maindata2$GDP)
maindata3<-maindata2%>%filter(Cons==1|Cons==2|Cons==5)
table(maindata3$ForT,maindata3$Cons)
table(maindata3$ForT)
maindata3$decision<-ifelse(maindata3$ForT==0&maindata3$Cons==5,1,0)
table(maindata3$decision)
table(maindata3$ForT,maindata3$Cons)
xtabs(~ForT+Cons+Clnt,data = maindata3)
table(maindata3$traded)
maindata3$logGDPT<-log(maindata3$GDPT)
maindata3$lndistance<-log(maindata3$capdist)
maindata3$logcinctrgt<-log(maindata3$cinctrgt)
maindata3$logcincclnt<-log(maindata3$cinccln)
summary(maindata3$YEAR)
table(YEAR)
maindata3$post911<-ifelse(maindata3$YEAR>2001,1,0)
table(maindata3$post911)
summary(cincratio)
summary(gdpratio)
maindata4<-maindata3%>%filter(maindata3$AgSt==2|maindata3$AgSt==3)
maindata4<-as.data.frame(maindata4)
maindata4$cincratio<-(maindata4$cinccln)/(maindata4$cinccln+maindata4$cinctrgt)
maindata4$gdpratio<-(maindata4$GDP)/(maindata4$GDP+maindata4$GDPT)
table(maindata4$ForT,maindata4$Cons)
table(maindata4$ForT,maindata4$Clnt,maindata4$Cons)
table(maindata4$decision)
table(maindata4$ForT)
table(maindata4$YEAR)
hist(maindata4$YEAR,xlab = "Year",ylab ="",main = "PMSC Throughout Years")
maindata4$postCold<-ifelse(maindata4$YEAR>1991&maindata4$YEAR<2001,1,0)
table(maindata4$defwithclnt)


# Checking for who use PMSCs


hist(maindata4$clientcode, main = "Histogram of Users")
table(maindata4$clientcode)
smalldata<-maindata4%>%group_by(clientcode)%>%filter(n()>=10)
table(smalldata$clientcode)
smallerdata<-smalldata%>%select(clientcode,Serv)
table(smallerdata$clientcode)
smallestdata<-smallerdata%>%filter(Serv<=7)
smallestdata<-smallestdata%>%group_by(clientcode)%>%mutate(use_count=n())
table(smallestdata$clientcode)
hist(smallestdata$use_count)


# Mapping


library(sf)
library(sp)
library(rnaturalearth)
library(rnaturalearthdata)
library(ggplot2)

smallerdata<-smallerdata%>%group_by(clientcode)%>%mutate(use_count=n())
world_map<-ne_countries(scale = "medium",returnclass = "sf")%>%filter(admin!="Antratica")
target_crs<-"+proj=moll"
world_map2<-st_as_sf(x=world_map)
world_moll<-world_map2%>%st_transform(crs = target_crs)
smallerdata$iso3c<-countrycode(smallerdata$clientcode,origin = "cown",destination = "iso3c")
filter1<-smallerdata$iso3c
world_moll2<-world_moll%>%dplyr::filter(iso_a3%in%filter1)
world_moll2%>%left_join(smallerdata,by=c("iso_a3"="iso3c"))%>%
  ggplot()+
  geom_sf(aes(fill=use_count))+
  scale_fill_viridis_c(trans="sqrt",labels=scales::percent_format(scale = 1),breaks=c(1:5)^2)+
  theme_bw()+
  theme(panel.background = element_rect(fill = "aliceblue"),plot.title = element_text(hjust = 0.5))+
  labs(title = "PMSC Use",caption = "Data Source: Petersohn et al. (2022)")


# Statistical Work for the Article
## Here we use the base model and the model 4 with relative information (in a more disintegrated approach) to analyze and publish results
### This reflects much more realistic approach to decision-making procedure that mirrors the real life style which is mostly based on relative information
#### As a result, the article finds this point of view more useful and less cumbersome for its analysis


## More check models sensitivity & specificity


attach(maindata4)
options(scipen = 999)
library(caret)
library(pROC)
library(ROCR)
library(dplyr)


partition1<-createDataPartition(maindata4$decision,times = 1,p=0.7,list = F)
trainset<-maindata4[partition1,]
testset<-maindata4[-partition1,]
basemodel1s<-glm(decision~as.factor(AgSt)+combat+mconsult+mtrain+mopsupport+mlogist+intellig+secur+traded,data = trainset,family = binomial(link = "logit"))
basemodelpred<-(predict(basemodel1s,newdata = testset,type = "response")>0.5)%>%as.numeric%>%as.factor
levels(basemodelpred)
levels(truval)
truval<-testset$decision
caret::confusionMatrix(basemodelpred,as.factor(truval),positive="1")
modeltargets<-glm(decision~as.factor(AgSt)+combat+mconsult+mtrain+mopsupport+mlogist+intellig+secur+traded+mounttrgt+capdist+atop_defense+logcinctrgt+polscrtrgt+logGDPT+post911,data = trainset,family = binomial(link = "logit"))
modeltargetspred<-(predict(modeltargets,newdata = testset,type = "response")>0.5)%>%as.numeric%>%as.factor
caret::confusionMatrix(modeltargetspred,as.factor(truval),positive="1")
modelrelatives<-glm(decision~as.factor(AgSt)+combat+mconsult+mtrain+mopsupport+mlogist+intellig+secur+traded+mounttrgt+capdist+defwithclnt+cincratio+abspol+gdpratio+post911,data = trainset,family = binomial(link = "logit"))
modelrelativespred<-(predict(modelrelatives,newdata = testset,type = "response")>0.5)%>%as.numeric%>%as.factor
caret::confusionMatrix(modelrelativespred,as.factor(truval),positive="1")


# ROC curves

basemodelpredict<-predict(basemodel1s,newdata = testset,type = "response")
roc(testset$decision,basemodelpredict)
plot(roc(testset$decision,basemodelpredict))
targetmodelpredict<-predict(modeltargets,newdata = testset,type = "response")
roc(testset$decision,targetmodelpredict)
plot(roc(testset$decision,targetmodelpredict))
relativemodelpredict<-predict(modelrelatives,newdata = testset,type = "response")
roc(testset$decision,relativemodelpredict)
plot(roc(testset$decision,relativemodelpredict))


## Models in the Article


# Base Model (information limited to the PMSC)

basemodel1<-glm(decision~as.factor(AgSt)+combat+mconsult+mtrain+mopsupport+mlogist+intellig+secur+traded,data = maindata4,family = binomial(link = "logit"))
summary(basemodel1)
summary(basemodel1,vcov=vcovCL(basemodel1,Clnt,type = "HC0"))
summary(basemodel1,vcov=vcovHC(basemodel1,type = "HC0"))
modelbasemargins<-marginaleffects(basemodel1,vcov = T,conf_level = .95,type = "response")
summary(margins(basemodel1,type = "response",vce = "bootstrap"))
basemodelmargins<-summary(modelbasemargins)


# Target Model (Information related to the Target + Time Related)

modeltarget<-glm(decision~as.factor(AgSt)+combat+mconsult+mtrain+mopsupport+mlogist+intellig+secur+traded+mounttrgt+capdist+atop_defense+logcinctrgt+polscrtrgt+logGDPT+post911,data = maindata4,family = binomial(link = "logit"))
summary(modeltarget)
summary(modeltarget,vcov=vcovCL(modeltarget,Clnt,type = "HC0"))
summary(modeltarget,vcov=vcovHC(modeltarget,type = "HC0"))
modeltargetmargins<-marginaleffects(modeltarget,vcov = T,conf_level = .95,type = "response")
summary(margins(modeltarget,type="response",vce="bootstrap"))
targetmodelmargins<-summary(modeltargetmargins)


# Relative Model (Information related to the relative comparison + Time Related)

modelrelative<-glm(decision~as.factor(AgSt)+combat+mconsult+mtrain+mopsupport+mlogist+intellig+secur+traded+mounttrgt+capdist+defwithclnt+cincratio+abspol+gdpratio+post911,data = maindata4,family = binomial(link = "logit"))
summary(modelrelative)
summary(modelrelative,vcov=vcovCL(modelrelative,Clnt,type = "HC0"))
summary(modelrelative,vcov=vcovHC(modelrelative,type = "HC0"))
modelrelativemargins<-marginaleffects(modelrelative,vcov = T,conf_level = .95,type = "response")
summary(modelrelativemargins)
summary(margins(modelrelative,type="response",vce="bootstrap"))
relativemodelmargins<-summary(modelrelativemargins)


# Stargazer of models

stargazer(basemodel1,modeltarget,modelrelative,model.names = T,style = "ajps",type = "text",align = T,digits = 5,out = "Model-Summary.txt",min.max = T,dep.var.labels = "Hiring PMSC by Foreign Government for Itself",covariate.labels = c("International vs Local","Combat Service","Military Consulting",
                                                                                                                                                                               "Military Training","Military Op.Support","Military Logistic","Intelligence",
                                                                                                                                                                                "Security","Publicly Traded","Mountainous Rage","Capital Distant",
                                                                                                                                                                               "Defense Alliance","CINC-Target(log)","Polity2-Target","GDP-Target","Defense with Client","CINC Ratio",
                                                                                                                                                                               "Polity2 Score Diff.","GDP Ratio","Post 9/11"))

# Visual Representation


library(visreg)
par(mfrow=c(2,2))


visreg(basemodel1,scale = "response","combat",xlab="Combat Service",ylab="P(hiring PMSC for the state)",main="Base-Model")
visreg(modeltarget,scale = "response","combat",xlab="Combat Service",ylab="P(hiring PMSC for the state)",main="Target-Model")
visreg(modelrelative,scale = "response","combat",xlab="Combat Service",ylab="P(hiring PMSC for the state)",main="Relative-Model")
visreg(basemodel1,scale = "response","mtrain",xlab="Military Training Service",ylab="P(hiring PMSC for the state)",main="Base-Model")
visreg(modeltarget,scale = "response","mtrain",xlab="Military Training Service",ylab="P(hiring PMSC for the state)",main="Target-Model")
visreg(modelrelative,scale = "response","mtrain",xlab="Military Training Service",ylab="P(hiring PMSC for the state)",main="Relative-Model")
visreg(basemodel1,scale = "response","mlogist",xlab="Military Logistics",ylab="P(hiring PMSC for the state)",main="Base-Model")
visreg(modeltarget,scale = "response","mlogist",xlab="Military Logistics",ylab="P(hiring PMSC for the state)",main="Target-Model")
visreg(modelrelative,scale = "response","mlogist",xlab="Military Logistics",ylab="P(hiring PMSC for the state)",main="Relative-Model")
visreg(basemodel1,scale = "response","secur",xlab="Security Services",ylab="P(hiring PMSC for the state)",main="Base-Model")
visreg(modeltarget,scale = "response","secur",xlab="Security Services",ylab="P(hiring PMSC for the state)",main="Target-Model")
visreg(modelrelative,scale = "response","secur",xlab="Security Services",ylab="P(hiring PMSC for the state)",main="Relative-Model")
visreg(modelrelative,scale = "response","defwithclnt",xlab="Defensive Agreement with the target",ylab="P(hiring PMSC for the state)",main="Relative-Model")
visreg(modeltarget,scale = "response","polscrtrgt",xlab="Polity 2 Score of the target",ylab="P(hiring PMSC for the state)",main="Target-Model")
visreg(modelrelative,scale = "response","abspol",xlab="Absolute Difference in Polity 2 Score",ylab="P(hiring PMSC for the state)",main="Relative-Model")
visreg(modeltarget,scale = "response","capdist",xlab="Capital Distance",ylab="P(hiring PMSC for the state)",main="Target-Model")
visreg(modelrelative,scale = "response","capdist",xlab="Capital Distance",ylab="P(hiring PMSC for the state)",main="Relative-Model")


# Plotting AME Results

library(ggplot2)

averagebase<-avg_slopes(basemodel1,conf_level = .95,vcov = T,type = "response")
ggplot(averagebase,aes(y=term,x=estimate,xmin=conf.low,xmax=conf.high))+
  geom_pointrange()+
  theme_bw()+
  theme(plot.title = element_text(hjust = 0.5))+
  geom_vline(xintercept = 0,linetype="longdash",col="red")+
  xlab("Average Marginal Effect")+
  ylab("Independent Variables")+
  labs(title = "AME for Base Model")

averagetarget<-avg_slopes(modeltarget,conf_level = .95,vcov = T,type = "response")
ggplot(averagetarget,aes(y=term,x=estimate,xmin=conf.low,xmax=conf.high))+
  geom_pointrange()+
  theme_bw()+
  theme(plot.title = element_text(hjust = 0.5))+
  geom_vline(xintercept = 0,linetype="longdash",col="red")+
  xlab("Average Marginal Effect")+
  ylab("Independent Variables")+
  labs(title = "AME for Target Model")
  

averagerelative<-avg_slopes(modelrelative,conf_level = .95,vcov = T,type = "response")
ggplot(averagerelative,aes(y=term,x=estimate,xmin=conf.low,xmax=conf.high))+
  geom_pointrange()+
  theme_bw()+
  theme(plot.title = element_text(hjust = 0.5))+
  geom_vline(xintercept = 0,linetype="longdash",col="red")+
  xlab("Average Marginal Effect")+
  ylab("Independent Variables")+
  labs(title = "AME for Relative Model")


## Supplementary Material: Including Marginal Average Effects and Marginal Discrete Effects Table


modellist<-list("Base-Model"=basemodelmargins,"Model-Target"=targetmodelmargins,"Model-Relative"=relativemodelmargins)
modelsummary::modelsummary(models = modellist,stars = T,fmt = "%.5f",output = "MarginalEffectsSupp.docx",title = "Table: Marginal Effects of Client State Decision")
stargazer(as.data.frame(maindata4),align = T,digits = 5,style = "ajps",out = "Descriptive Table for Supplementary Material.txt",type = "text",min.max = T)
sessionInfo()%>%capture.output(file = "session_info.txt")

